library(tidyverse)
library(readxl)
library(ggplot2)
library(zoo) # for function rollapply
library(scales) # for rescaling btw 0 and 1
library(plotly) # for interactive plot https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html
library(hrbrthemes) # for https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html
library(htmlwidgets) # for https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html
library(viridis) # for https://www.r-graph-gallery.com/stacked-area-chart-plotly.html
library(gganimate) # for animated plot https://www.r-graph-gallery.com/287-smooth-animation-with-tweenr.html
library(leaflet) # for map
library(magrittr)
library(htmlwidgets)
library(rgdal)
#number of fishers
fishers<-read.csv(file.path("data/number_age_parttime_fulltime_fishers.csv"), sep=",")
#study municipalities list
municip<-read.csv(file.path("C:/Users/sigrid.engen/github/nor-prep/prep/administrative/komlist.csv"), sep=";")
#population in study municipalities
population<-read.csv(file.path("C:/Users/sigrid.engen/github/nor-prep/prep/administrative/data/population_size_study_municipalities_1986_2018.csv"),sep="," )
#spatial layer of municipalities in Norway
spat_municip<-readOGR(dsn="C:/Users/sigrid.engen/github/nor-prep/prep/administrative/raw/shapefiles/NO_AdminOmrader", "NO_AdminOmrader_pol")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\sigrid.engen\github\nor-prep\prep\administrative\raw\shapefiles\NO_AdminOmrader", layer: "NO_AdminOmrader_pol"
## with 428 features
## It has 7 fields
fishers_0<-left_join(fishers, municip, by = c("municip_numb" = "Komnum")) %>%
dplyr::select(-c(X, Kid, municip_name)) %>%
rename(municip_name = "Name")
# summarise fishers across the age variable
fishers_1<- fishers_0 %>%
group_by(municip_numb,
municip_name,
year) %>%
summarise(main_occup=sum(main_occupation)) %>%
ungroup()
## `summarise()` regrouping output by 'municip_numb', 'municip_name' (override with `.groups` argument)
# calculate moving average
### the rollapply function places the calculated value on the row = year three.
### We want it to be placed year 6 - we want the moving average of the 5 previous years to
### compare this with the value the 6th year
fishers_2 <- fishers_1 %>%
group_by(municip_numb) %>%
mutate(five_yr_roll_aver = rollapply(main_occup, 5, mean, fill=NA)) %>%
dplyr::select(-main_occup) %>%
mutate(year = year + 3) %>% #recoding year variable - have to remove count and add back in in next step so the order is correct
filter(!year %in% c("2020", "2021", "2022")) %>%
ungroup()
# add values back to ensure that order of year and values are correct and calculate scores
fishers_3<- left_join(fishers_1, fishers_2, by = c("municip_numb", "municip_name", "year")) %>%
filter(!is.na(five_yr_roll_aver)) %>% # remove NAs
mutate(main_occup_score = main_occup/five_yr_roll_aver) %>% #divide value the 6th year with average the 5 preceding years
mutate(main_occup_score_1 = replace(main_occup_score, main_occup_score >= 1, 1)) # scores that are larger than 1 because the number of fishers are larger than the average the 5 preceeding years are given a value of 1
ggplot(fishers_3)+
geom_col(aes(x = year, y = main_occup_score_1), fill = "dodgerblue4")+
facet_wrap(municip_name ~ .) +
labs(y = "Sustainability score", x = "") +
ggtitle("")+
theme(legend.position = "none",
axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
plot.title = element_text(size=3, hjust=0.5, vjust=0.5),
axis.ticks.x = element_line(size = 0.5),
axis.ticks.length = unit(0.2, "cm"),
text=element_text(family="sans"))
fishers_1_p<- fishers_0 %>%
group_by(municip_numb,
municip_name,
year) %>%
summarise(sum(part_time_occupation)) %>%
dplyr::rename(part_occup = 'sum(part_time_occupation)') %>%
ungroup()
## `summarise()` regrouping output by 'municip_numb', 'municip_name' (override with `.groups` argument)
fishers_2_p <- fishers_1_p %>%
group_by(municip_numb) %>%
mutate(five_yr_roll_aver = rollapply(part_occup, 5, mean, fill=NA)) %>%
dplyr::select(-part_occup) %>%
mutate(year = year + 3) %>%
filter(!year %in% c("2020", "2021", "2022")) %>%
ungroup()
fishers_3_p<- left_join(fishers_1_p, fishers_2_p, by = c("municip_numb", "municip_name", "year")) %>%
filter(!is.na(five_yr_roll_aver)) %>%
mutate(part_time_occup_score = part_occup/five_yr_roll_aver) %>%
mutate(part_time_occup_score_1 = replace(part_time_occup_score, part_time_occup_score >= 1 & part_time_occup_score < Inf, 1))
ggplot(fishers_3_p)+
geom_col(aes(x = year, y = part_time_occup_score_1), fill = "dodgerblue4")+
facet_wrap(municip_name ~ .) +
labs(y = "Number of fishers", x = "") +
ggtitle("test")+
theme(legend.position = "none",
axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
plot.title = element_text(hjust=0.5, vjust=0.5),
axis.ticks.x = element_line(size = 0.5),
axis.ticks.length = unit(0.2, "cm"))
## Warning: Removed 5 rows containing missing values (position_stack).
fishers_all<- fishers_0 %>%
group_by(municip_numb,
municip_name,
year) %>%
summarise_at(c("main_occupation", "part_time_occupation"), sum) %>%
mutate(total_fishers = main_occupation+part_time_occupation) %>%
ungroup()
fishers_all_2 <- fishers_all %>%
group_by(municip_numb) %>%
mutate(five_yr_roll_aver = rollapply(total_fishers, 5, mean, fill=NA)) %>%
select(-c(main_occupation:total_fishers)) %>%
mutate(year = year + 3) %>%
filter(!year %in% c("2020", "2021", "2022")) %>%
ungroup()
fishers_all_3<- left_join(fishers_all, fishers_all_2, by = c("municip_numb", "municip_name", "year")) %>%
filter(!is.na(five_yr_roll_aver)) %>%
mutate(occup_score = total_fishers/five_yr_roll_aver) %>%
mutate(occup_score_1 = replace(occup_score, occup_score >= 1, 1))
ggplot(fishers_all_3)+
geom_col(aes(x = year, y = total_fishers), fill = "dodgerblue4")+
facet_wrap(municip_name ~ .) +
labs(y = "Number of fishers", x = "") +
ggtitle("")+
theme(legend.position = "none",
axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
plot.title = element_text(hjust=0.5, vjust=0.5),
axis.ticks.x = element_line(size = 0.5),
axis.ticks.length = unit(0.2, "cm"))
### Plot total number of fishers by municipality (without Tromsø)
fishers_all_3_noTromso <- fishers_all_3 %>%
filter(!municip_name %in% c("Tromso"))
ggplot(fishers_all_3_noTromso)+
geom_col(aes(x = year, y = total_fishers), fill = "dodgerblue4")+
facet_wrap(municip_name ~ .) +
labs(y = "Number of fishers", x = "") +
ggtitle("")+
theme(legend.position = "none",
axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
plot.title = element_text(hjust=0.5, vjust=0.5),
axis.ticks.x = element_line(size = 0.5),
axis.ticks.length = unit(0.2, "cm"))
ggplot(fishers_all_3)+
geom_col(aes(x = year, y = occup_score_1), fill = "dodgerblue4")+
facet_wrap(municip_name ~ .) +
labs(y = "Sustainability score", x = "") +
ggtitle("")+
theme(legend.position = "none",
axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
plot.title = element_text(hjust=0.5, vjust=0.5),
axis.ticks.x = element_line(size = 0.5),
axis.ticks.length = unit(0.2, "cm"))
### Plot distribution of all scores
ggplot(fishers_all_3, aes(occup_score_1)) +
geom_bar(fill = "blue") +
scale_x_binned()
### Plot distribution of scores by municipality
ggplot(fishers_all_3, aes(occup_score_1)) +
facet_wrap(municip_name ~ .) +
geom_bar(fill = "blue") +
scale_x_binned()+
theme(axis.text.x = element_text(size = 5))
## Plot total number of fishers - interactive plot
fishers_northernNorw<- fishers %>%
group_by(year) %>%
summarise_at(c("main_occupation", "part_time_occupation"), sum) %>%
mutate(total_fishers = main_occupation+part_time_occupation) %>%
ungroup()
# Usual area chart
p <- fishers_northernNorw %>%
ggplot(aes(x=year, y=total_fishers)) +
geom_area(fill="#69b3a2", alpha=0.5) +
geom_line(color="#69b3a2") +
ggtitle("Utvikling i antall fiskere i Nord-Norge 1983-2018, Data: Fiskeridirektoratet") +
ylab("Antall fiskere") +
xlab("Ã…r") +
theme_ipsum() +
theme(plot.title = element_text(family = "arial", size=12, face= "italic"),
axis.title.x = element_text(family = "arial", size=10),
axis.title.y = element_text(family = "arial", size=10),
axis.line = element_line(size = 2, colour = "grey80"))
# Turn it interactive with ggplotly
p <- ggplotly(p)
p
# save the widget
#saveWidget(p, "figs/ggplotlyAreachart_Total_Fishers_NN.html", selfcontained = TRUE)
fishers_northernNorw_long<- fishers_northernNorw %>%
select(-c(total_fishers)) %>%
pivot_longer(cols = c("main_occupation", "part_time_occupation"),
names_to = "fishers",
values_to = "n")
# Plot
p_2 <- fishers_northernNorw_long %>%
ggplot(aes(x=year, y=n, group=fishers, color=fishers)) +
geom_line() +
geom_point() +
scale_color_viridis(discrete = TRUE) +
ggtitle("Utvikling i antall fiskere i Nord-Norge") +
theme_ipsum() +
ylab("Antall fiskere") +
xlab("Ã…r")+
transition_reveal(year) +
theme_set(theme_gray(base_size = 20, base_family = 'sans' )) +
scale_color_manual(values=c('#008B8B','#000000'))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
animate(p_2, duration = 5, fps = 20, width = 500, height = 500, renderer = gifski_renderer())
#anim_save("output.gif")
# Save at gif:
#anim_save(p_2, "figs/287-smooth-animation-with-tweenr.gif")
#getwd()
fishers_age_NN<- fishers %>%
group_by(year,
age) %>%
summarise_at(c("main_occupation", "part_time_occupation"), sum) %>%
mutate(total_fishers = main_occupation+part_time_occupation) %>%
ungroup()
#calculate age distribution in percentages
fishers_age_NN_percentages<- fishers_age_NN %>%
group_by(year) %>%
mutate(percentage_main = main_occupation/sum(main_occupation)*100) %>%
mutate(percentage_part = part_time_occupation/sum(part_time_occupation)*100) %>%
ungroup()
# Plot age distribution main occupation
p <- fishers_age_NN_percentages %>%
ggplot( aes(x=year, y=percentage_main, fill=age, text=age)) +
geom_area( ) +
scale_fill_viridis(discrete = TRUE) +
theme(legend.position="none") +
ggtitle("Aldersfordeling (%) blant fiskere med fiske som hovedyrke, Data: Fiskeridirektoratet") +
theme_ipsum() +
theme(legend.position="right")+
theme(plot.title = element_text(family = "arial", size=12, face= "italic"),
axis.title.x = element_text(family = "arial", size=10),
axis.title.y = element_text(family = "arial", size=10),
axis.line = element_line(size = 2, colour = "grey80"))
# Turn it interactive
p <- ggplotly(p, tooltip = c("text", "x", "fill"))
p
# save the widget
# library(htmlwidgets)
# saveWidget(p, file=paste0( getwd(), "/HtmlWidget/ggplotlyStackedareachart.html"))
# Plot age distribution main occupation
p <- fishers_age_NN_percentages %>%
ggplot(aes(x=year, y=percentage_part, fill=age, text=age)) +
geom_area( ) +
scale_fill_viridis(discrete = TRUE) +
theme(legend.position="none") +
ggtitle("Aldersfordeling (%) blant fiskere med fiske som biyrke, Data: Fiskeridirektoratet") +
theme_ipsum() +
theme(legend.position="right")+
theme(plot.title = element_text(family = "arial", size=12, face= "italic"),
axis.title.x = element_text(family = "arial", size=10),
axis.title.y = element_text(family = "arial", size=10),
axis.line = element_line(size = 2, colour = "grey80"))
# Turn it interactive
p <- ggplotly(p, tooltip = c("text", "x", "fill"))
p
# save the widget
# library(htmlwidgets)
# saveWidget(p, file=paste0( getwd(), "/HtmlWidget/ggplotlyStackedareachart.html"))
##################################
### Prepare attribute data #######
##################################
# Check correspondence btw the fishers file and study municip file in terms of municipal numbers
e<-anti_join(fishers, municip, by =c("municip_numb" = "Komnum")) # no municip numbers in fishers not in municip
f<-anti_join(municip, fishers, by =c("Komnum" = "municip_numb")) # noncoast is the only variable in municip not in fishers
# extract fishers data for 2018 for plotting on the map
fishers_northernNorw_2018<- fishers %>%
filter(year =="2018") %>%
group_by(year,
municip_numb,
municip_name) %>%
summarise_at(c("main_occupation", "part_time_occupation"), sum) %>%
mutate(total_fishers = main_occupation+part_time_occupation) %>%
ungroup()
##################################
### Prepare spatial data #########
##################################
#Change municipal number of Harstad from 1901 to the correct 1903 in the spatial file over municipalities in NOrway. THe muncipality was merged with Bjarkjøy in 2013 and went from number 1901 to 1903. The polygon in the shapefile reflects this merger (Harstad includes Bjarkøy) but the number is wrong.
spat_municip@data[spat_municip@data$KOMM == "1901", "KOMM"] <- "1903"
#extract subset of data (studymunicipalities from komlist.csv) from spatialpolygonsdataframe
muni_sub <- spat_municip[spat_municip$KOMM %in% c(1804L, 1805L, 1811L, 1812L, 1813L, 1815L, 1816L, 1818L, 1820L,
1822L, 1824L, 1827L, 1828L, 1832L, 1833L, 1834L, 1835L, 1836L,
1837L, 1838L, 1840L, 1841L, 1845L, 1848L, 1849L, 1850L, 1851L,
1852L, 1853L, 1854L, 1856L, 1857L, 1859L, 1860L, 1865L, 1866L,
1867L, 1868L, 1870L, 1871L, 1874L, 1902L, 1903L, 1911L, 1913L,
1917L, 1919L, 1920L, 1923L, 1924L, 1925L, 1926L, 1927L, 1928L,
1929L, 1931L, 1933L, 1936L, 1938L, 1939L, 1940L, 1941L, 1942L,
1943L, 2002L, 2003L, 2004L, 2012L, 2014L, 2015L, 2017L, 2018L,
2019L, 2020L, 2022L, 2023L, 2024L, 2025L, 2027L, 2028L, 2030L
), ]
#add number of fishers to spatial layer of study municipalities
muni_fishers<- sp::merge(muni_sub, fishers_northernNorw_2018, by.x = "KOMM", by.y= "municip_numb")
#add names of municipalities from list to spatial layer of study municipalities
muni_fishers_1<- sp::merge(muni_fishers, municip, by.x = "KOMM", by.y= "Komnum")
#transform projection of this spatial layer to fit leaflet and open street map
#https://stackoverflow.com/questions/34769845/how-to-plot-polygons-on-leaflet-using-r
muni_fisher_transform <- spTransform(muni_fishers_1, CRS("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs"))
##################################
###Build map #####################
##################################
#https://www.theanalyticslab.nl/polygon-plotting-in-r/
#https://rstudio.github.io/leaflet/colors.html
#https://www.r-graph-gallery.com/183-choropleth-map-with-leaflet.html
#Define cut points for the colorbins refering to number of fishers
cuts <- c(0.0, 1, 50, 100, 150, 200, 250, 300, 350, 400)
# Choose a color palette and assign it to the values
colorbins <- colorBin("Blues", domain = muni_fisher_transform$total_fishers, bins = cuts)
# Prepare the text for tooltips:
mytext <- paste(
"Kommune: ", muni_fisher_transform@data$Name,"<br/>",
"Antall fiskere: ", muni_fisher_transform@data$total_fishers) %>%
lapply(htmltools::HTML)
# Show the number of fishers by municipality on a Leaflet map
map<-leaflet(muni_fisher_transform) %>%
setView(lng = 20, lat = 68, zoom = 5) %>%
addTiles() %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addPolygons(stroke = TRUE,
color = "white",
weight="1",
smoothFactor = 0.3,
fillOpacity = 0.7,
fillColor = ~colorbins(muni_fisher_transform$total_fishers),
label = mytext,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "13px",
direction = "auto"
))
# Add a legend
map_with_legend <- map %>% addLegend(pal = colorbins,
values = muni_fisher_transform$total_fishers,
opacity = 0.7, title = "Number of fishers 2018", position = "bottomright")
map_with_legend
###########################
## SAVE OUTPUT HTML FILE ##
###########################
# Save output as HTML widget (or incorporate into Shiny / Flexdashboard)
#saveWidget(map_with_tooltip, file="Elderly residents in Utrecht.html")